# Estimate P-value threshold
# 2013-10-29

# Rscript FDR_filter.R [target aso file with Pval] [random aso file with Pval] [filtered output file]

# False Discovery Rate threshold
FDR_level <- 0.05

targetFile <- commandArgs(trailingOnly =T)[1]
randomFile <- commandArgs(trailingOnly =T)[2]
outputFile <- commandArgs(trailingOnly =T)[3]

target <- read.table(targetFile, header=F, sep="\t")
random <- read.table(randomFile, header=F, sep="\t")
targetP <- target$V6
randomP <- random$V6

randomP <- sort(randomP)
Threshold <- randomP[as.integer(0.05*length(randomP))]

#mixp <- c(targetP, randomP)
#mixp_uniqu <- unique(mixp[mixp < 0.05 & mixp > 0])
#FDRs <- sapply(mixp_uniqu, function(x) sum(randomP < x) / sum(mixp < x))
#Threshold <- max(FDRs[FDRs <= FDR_level], na.rm=T)

totalNum <- length(targetP)
PvalCandidate <- sum(targetP < 0.05, na.rm=T)
cat(paste(as.character(PvalCandidate), "/", as.character(totalNum), "\n"))

cat(paste("Threshold:", as.character(round(Threshold,5)), "\n"))
totalSig <- sum(targetP < Threshold, na.rm=T)

cat(paste(as.character(totalSig), "/", as.character(PvalCandidate), "\n"))

write.table(subset(target, V6 < Threshold), file=outputFile, quote=F, sep="\t", row.names=F, col.names=F)

